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SUMMARY 


A numerical study of finite-strength, isentropic pressure waves trans- 
verse to the axis of a circular cylinder is made for the radial resonant 
mode. The waves occur in a gas otherwise at rest, filling the cylinder. 

A method of characteristics is used for the numerical solution. 

For small but finite amplitudes, calculations indicate the existence 
of waves of permanent potential form. For larger amplitudes, a shock is 
indicated to occur. The critical value of the initial amplitude parameter 
in the power series is found to be 0.06 to 0.08, under various types of 
initial conditions. 


INTRODUCTION 


Combustion instability involving large pressure oscillations ' in the com- 
bustion chambers of rocket and jet engines has been observed. In (1), (2) and 
(3) , these oscillations were analyzed on the assumption that oscillations were 
so weak that the classical acoustic theory is valid. However, a propagating 
plane wave of finite strength will eventually develop into a shock. Thus, 
one might expect that the high observed pressures imply the presence of shocks. 
Actually, in a cylindrical case, a shock may possibly fail to appear owing to 
the wave-scattering effect of the chamber wall. Moore and Maslen [4] studied 
finite transverse waves in a circular cylinder, assuming that the disturbance 
velocity potential could be expanded in a power series in amplitude. They 
found solutions for the first three approximations to be periodic in time, 
whereas only the first approximation is purely periodic in the plane-wave case. 
They concluded that finite-strength, shock-free oscillations in a cylinder 
were possible for some range of amplitude. 

The present study employs numerical solutions of the dynamic equations to 
find a critical value of the amplitude, up to which the power series converges 
for the radial cylindrical mode, and there exists a permanent strong trans- 
verse oscillation without a shock. 
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FORMULATION OF THE PROBLEM 


Differential Equations 

An isentropic, irrotational fluid motion is governed by the following 
equations : 


Continuity 

p fc + V • (p_g) = 

0 

(la) 


_ 2 

-2 


Momentum 

<p + -V + — T £ = 
t 2 y-1 p 

a 

Y-1 

(lb) 

State isentropy 

Y 

P 02 P 


(lc) 

Velocity potential 

a = v< f) 


(Id) 


where a is the speed of sound in the absence of any oscillation. By elimin- 
ating p, p and from eqs . (1), the following equation for velocity potential 
is obtained: 


it - (52 - + I ♦? + 2 ^ <r V + £ »ee> 

r 


+ { 97 + ¥ ^r 2^2 = 0 


r r^ 


( 2 ) 


In this report, we shall consider only axisymmetric cases (-r^- = 0) . 
Also, for simplicity, y - 1 ~ 0 is assumed. Thus, equation (2) becomes 


-2 

d> — (a 2 — d> 2 )<j> + 2d) <b - — — d> = 0 

T tt v f r^rr Vrt r *r 


(3) 


Denoting dimensional variables with primes, nondimensional variables are 
defined : 


r = 


r 

R 


t = t’a/R 


d>' 

d> = i— 
v aR 
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where R is the radius of the cylinder. The resulting non-dimensional equa- 
tion is 


i - (1 - 4>2) <j> +2<f>cj> - — <f> = 0 

tt Y r T rr Y r Y rt r Y r 


(4) 


The following boundary conditions specify no flow through the cylinder wall 
and no flow at the cylinder axis : 


4> r (1 , t ) = 4> r (0,t) = 0 


(5) 


The corresponding dimensionless pressure is: 

4> 2 

’(+t + HT> . 

p = e - 1 


( 6 ) 


where p = p'/(p a^). p is the density of the fluid in the chamber when there 
o o 

is no disturbance, and p is zero if there is no disturbance in the chamber. 


NUMERICAL METHOD 


Characteristics 


Equation (4) proves to be a hyperbolic equation. Therefore, a method of 
characteristics can be used for numerical solution (Ref. 5.) . 


Let 


u = d> , v = d> 

y t ’ y r 


(7) 


Substituting equations (7) into equation (4) yields: 


u -(1 - v^)u + 2v u = — v. 
t r r r 


( 8 ) 


Directly from equation (7), 


u = v 
r t 


(9) 


If we now introduce new characteristic coordinates n , then the first 
derivatives with respect to r and t become: 


3 


_3_ 

3r 


(10a) 


r 3 , 3 

35 + n r 3n 


_3_ _ _ _3 . _3_ 

3t H 35 n t 3n 


and equations (8) and (9) become 


(10b) 


5 t u 5 + Vn - (1 - v2) (S r v 5 + Vn ) + 2v(5 r u e + \ u n> ' r v 


5 r u 5 + Vn ’ ( V{ + "tV * 0 


Thus 


n + 2vn , -d - V 2 )n 


u 


u (5 + 2v5 )-(l - v 2 ) 5 V - — V 

t r r 


n 

_ 

5 t r r 5 r 

n , - n 


V 


5 u_ - 5 v_ 

r t 


n 


L r 5 t 5 


In order to make u and v indefinite, we have to have 

n n 


n 2 (l-v 2 ) - n t (n t + 2vn r ) = 0 (11) 

and 


n t [(? t + 2v? r ^ U 5 " ( 1_v2 ^ r v ^ _ r v ]"(l-v 2 )n r (? r u^ - 5 t v = 0 (12) 


But 


11 1 _ _ 3r . 

n “ 3t ; 

r n 


and thus equation (11) becomes: 


f ) - 2v f ) - (1-v 2 ) - 0 

n n 
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which yields the result 


|f) = -(1-v) or (1 + v) (13) ,(14) 

n 

The corresponding form of equation (12) is 


If) [U t + 2vC r )u ? - (l-v 2 K r v 5 - | v] + (1-v 2 ) (5^ - C t v 5 ) = 0 (15) 
If -|f-) = -(1-v) ^ 0, eq. (15) becomes 

O t 

n 

(s t + 2v^ r )u^ -(l-v 2 )5 r v^ - | V -(1 + v)(C r U ? - C t v^) = o 
Thus, [£ t - (1 - v) C r ] [u^ + (1 + v)v^] = | v 
And since 


li) _ i£) + i£) K) „ 1£) _ (1 _• v ) r 

3t ; 9t J 9t ; 9r ; 9t ; U v; Sr 

n r n t r 


we have the characteristic equation 


r 9us \ 3v> 1 

8t> l 3? } + (1 + v) 3? 1 " r v 

n n n 


or 



+ (1 + v) |f) 


1 

r 


v 


n n 

9 r 

If, on the other hand, — ) = 1 + v ^ 0, then 

O t 

c [u - (1 - v)v ] + £ [2vu - (1-v ) v + (1 - v)u ] = I V 

t K t, x E, E, £ r 

and in a similar manner, we obtain the second of the pair of characteristic 
equations : 
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|^+(l+v)|^ = iv on 1^- = - (1 - v) i 0 (16a) 

- (1 - v) ^ - v on = (1 + v) (16b) 

3t dt r ot 

The pressure field is, (from Eq. (6)) 

p = exp[-(u + 7p)] ~ 1* (17) 


Numerical Calculation (Ref. 6) 



Fig. 1. Finding a new mesh point. 
Characteristics - AD and - BD are 
of opposite families. 


In all cases, the differential 
coefficients are replaced by the 
divided differences over an interval, 
and other quantities are replaced by 
the arithmetic mean over the inter- 
val. For example, over the range 
AD in Fig. 1, the equation 
dr 

— = 1 + v will be replaced by the 
equation r Q -r A = [1 + - t A ) 

where [1 + v]^ signifies an arithme- 
tic mean over the interval A and D 
i.e. , 1 + 1/2 (v + v ) . 


Suppose we know all values v, u, t and r at points A and B; then we can 
find these values at point D as follows: Eq. (16) is replaced by 


r D ‘ r A = [1 + V] AD (t D ' t A ) 


r D r B t 1 V ] BD (t D tg) 


u D - n A + [v - Hajj <v d - r A ) = (t D - t A ) fl i vdt* 
U D “ U B + *- V + 1 ^BD ^ V D ” V B^ " ^D ~ t B' > ^B 7 Vdt 


And from eqs . (18) and (19), 


For the calculation of this integral, see Appendix A. 


(18) 

(19) 

( 20 ) 
( 21 ) 
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( 22 ) 


C D " 


r B r A + + 1 '*AD ~ t B^ V ~ 1 ^BD 

+ 1 ^AD “ t v " H BD 


r D = r A + t v + 1]^^ - t A > 


Also, from eqs . (20) and (21), 


(23) 


V D * 


/“ i vdt x (t -t ) - (t -t )/“ | vdt + V [v-1] -v B [v+1] BD + u -u 
A B 


[v-1] 


AD 


[v+1] 


(24) 


BD 


U D = U A " [V “ 1] AD (V D ' V + r Vdt X (t D ’ t A ) 

A 


(25) 


Thus, we can find the values of u, v, r and t at the new point. Furthermore, 
using the foregoing equations, if we know all the information at the points 
1, 2, ........ , 10 on line A-A' of Fig. 2, we can calculate new values at 

the points 1', 2, , 10" on line BB'. 



Fig. 2. Grid of characteristics illustrating the process of calculation from 
the known state AA' to the new state BB' . 
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Shock Occurrence 



Fig. 3. A shock occurrence, where two characteristics of the same family 

cross. (At point P, 2 characteristics of family (I) are crossing.) 


A shock is taken to be indicated by the intersection of two characteris- 
tic lines of the same family. Following such an occurrence, the development of 
a shock will obey the Rankine-Hugoniot conditions. In this paper, the 
object is to find the possibility of the shock formation as a function of 
wave amplitude. Therefore, the calculations do not consider the propagation 
of the shock after it has first been indicated by coalescence of character- 
istics. However, it may be noted that, in our method, the time level at each 
stage of calculations is not the same. Therefore, in order to describe the 
wave shape at the time the shock is first indicated, additional calculations 
are needed to fill the empty (t,r) space up to that time. Referring to Fig. 

4, the upper right hand side (cross hatched part) must be filled in. 



Fig. 4. Distribution of grid points when shock is indicated. 
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Initial Conditions 


Equation (4) is non-linear; therefore, the choice of initial condition 
is very important for this particular study. For example, if a function 
with singularities in the first derivatives is picked as an initial velocity 
distribution, then a shock may occur very early at a very small magnitude. 

In order to provide a suitable variety of initial conditions, the following 
three problems are posed: 

(a) Use the power-series solution (see ref. [4] and also Appendix A) 
to the 2nd order approximation as an initial condition 


<t> = e<p 


( 1 ) 


+ £“ 


( 2 ) 


(26) 


and perform calculations for various values of initial amplitude e. 

(b) Use the solution of the power series to the 2nd order approximation 
with the velocity magnitude = 0.1 as initial condition and put the velocity 
disturbance at the outer boundary as a function of pressure at the wall, 
v wall = ~ °P ( ref • [3]) : 


v wall = " 0 * P (27) 

This gives a gradual amplification of the transverse oscillation and the rate 
of this magnification depends on the coefficient a. This case corresponds 
to the response of a vigorous burning zone near the wall. This response may 
be considered to provide a simple model of combustion instability. 

(c) Start with the procedure of case (b) , but after the magnitude of 
the oscillation reaches some prescribed value, then stop the amplification 
(set a=0) and watch the oscillation proceed. By this technique, the direct 
effect of the wall amplification condition on the shock formation may be 
minimized. Also, such an amplification cut-off would represent the limiting 
effects of unspecified damping mechanisms. Further, using small values of a 
should also minimize any direct effects of amplification rates on the waves.. 


Error Estimation 

There are no constraints imposed on the spatial mesh size, and A t is 
automatically determined once the spatial grid size is given. 

By using a Taylor series analysis of the equations, truncation errors 
for each equation can be estimated. 

If) - v-l (13) 

n 
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The Taylor series is 


3r, 


1 3 2 r, 


r i + i = r i + ^ . At+ i^> . At ' + o*t 3 ) 

n,i n,i 


3r. a , 1 3 2 r s 

r i ~ r i+l 3t . 2 3t2 ) 

n , i+l n » i+l 


At 2 + 0(At 3 ) 


3 2 r 


1 r 3r\ . 3r. , . ^ , 1 r 3 2 r. 

r i+l ~ r i = 2 3t . + 3t } . J At + 4*3t2* . 3 t 2 

n , i n , i+l n , i n , i+l 


) ]a t 2 + 


Applying Eq. (13) for both points i and i+l, we get 


3rs _ , 

at 5 . v i+i 1 
n,i+l 


3r. 

i? . ' v i - 1 

n,i 


We add above two equations and divide by 2: 


5tft> , +1 + f> .1- I''’ 1 '. , + , 

n , i+l n , i i , i+l 

Comparing Eq. (28) and (29) gives the truncation error: 


TE = I r 

4 L 3t2 ; 


n,i 


a£r 

3t 2 


) ]At + 0 (A t 2 ) 
rii+1 


Integrating eq. (16a) gives 


u 1 -u. + [-rv + v] - [-r v + v] = /. — vdt 

i+l l 2 . . - 2 . lr 

i+l l 


U i+1 


U i + (u i+l - V i H I (v i+l + V + 1] = f i +1 r V dt 


whence : 


U i+1 


u . + (v - v. ) [1 + v] . = / d+d — v dt 

l i+l l l , i+l l r 


(28) 
0 (A t 3 ) 


(29) 

(16a) 
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Truncation error for this equation is only that for the integration. (See 
Appendix B) 

For the integration, the assumption is made that the velocity profile 
between 2 grids (i & i + 1) is linear, i.e.. 


v 

o 


v. + 
1 



V . 
1 


Ar 


1 


But a Taylor Series expansion yields 



Fig. 5. A characteristic line 
between 2 grid points i and i+1. 
0 is a point on this character- 
istic line . 


v ... -v. 

, i+l i . 

v = v. H : Ar, 

o l Ar 1 


1 3 z v. , . 

2 '^2 )Ar l (Ar ^ r l ) 


+ 0 (A r 3 ) . 


So the error in the velocity due to the 
linear assumption is 


1 pi 2,. 

E ^ 0[j(^|) * Ar x (Ar - Ar )] 
o 


nr l,8 2 v. 
< ° 8 9r2 


Ar 2 ] 


max 


where (f^f) 


is the maximum value of 


max 


on the characteristic line between 


8 2 v 
3r2 

i and i+1. As is shown in Appendix A, 
the integration was carried out as 
follows . 


+ +1 “ v dt = / 1+1 i ? i!i rr . T +_ T / i+1 2. dr 

1 r i r [v+l] 11+1 [v+l] il+1 r r 

Thus, the resulting truncation error is of order: 


T.E. 


orl (li) _J 

Ul 8 t '3r2 ; [v+1] 

max 


Ar 3 , 


i , i+1 
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RESULTS 


For all calculations, 51 grid points are chosen. Thus, A r is 0.02 
initially. 

Case (a) . - The following calculations were performed for initial 
conditions given by two terms of a power series in e (eq. (26)): 


Table 1: Shock occurrences for various initial wave magnitudes. 


Initial A 
Magnitude 

e 

Shock 

At Time 
(K/I o ) 

0.1 

(4.485 x 10“ 2 ) 

No 

— 

0.175 

(7.849 x 10" 2 ) 

No 

— 

0.195 

(8.746 x 10“ 2 ) 

Yes 

13.5 

0.200 

(8.970 x 10“ 2 ) 

Yes 

5.2 

0.250 

(1.1213 x 10 -1 ) 

Yes 

2.2 

0.300 

(1.3456 x 10 -1 ) 

Yes 

1.5 

0.350 

(1.5698 x 10 _1 ) 

Yes 

1.3 


Fig. 6 shows the behavior of a wave with initial magnitude 0.30 (e =0.1346). 
From this figure and with the table above, it is clear that as time goes on, 
a sharp wave front appears for the waves which 'have comparatively large 
initial magnitudes and the crests grow, eventually being followed by shock 
formation. 

Changes of wave magnitude with time for each initial magnitude is descri- 
bed in Fig. 7. By magnitude is meant the highest wave amplitude during a 
half-cycle of wave oscillation. Every wave increases its magnitude at first, 
but lower magnitude waves often reach a certain amplitude, decay almost to 
their own initial conditions. Waves repeat this process with a period almost 
equal to three times that of the wave oscillation period. 

The amplitude of this magnitude change is of order (eg) 3 , i.e., 3rd order 
in the power series for velocity. Since there is no damping factor like 
viscosity in this problem, the difference between initial condition and the 


Definition of Magnitude (see also Appendix B): Magnitude = 0.581864 (eg), 

where g is the first zero of first order Bessel function, i.e., J^(g)=0. 
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exact solution for the differential equation does not attenuate and presumably 
causes these fluctuations of wave magnitude. Fig. 8 indicates the possibility 

Table 2: Order of the change of wave magnitude. (eB) 3 is the order of the 

3rd order approximation of the power series solution for velocity. 


Initial Magnitude 

(eB) 3 

A (Magnitude of Wave) 

0.1 

0.0051 

~ 0.0025 

0.175 

0.0272 

~ 0.025 

0.195 

0.0378 

~ 0.045 . 

0.200 

0.0407 

~ 0.050 


of the existance of a wave magnitude above which the wave grows to a shock, 
but below which it oscillates in the cylinder forever. An initial wave magni- 
tude around 0.180 seems to be the borderline for this particular initial 
condition. (e value of the series is 0.08.) 

Due to the irregularity of this oscillation, we failed in the attempt to 
find a repetition of identical wave patterns after a certain number of cycles 

for an initial magnitude of 0.175. However, as seen in Fig. 9, the wave shape 

of magnitude 0.175 after 16 cycles still is not steep at all. At least this 
helps to justify our conclusion that no shock occurs. 

Convergence and stability of the finite difference scheme were also 
checked to determine if round-off or truncation error could, be the source of 
fluctuations of wave magnitude. In Fig. 10, comparison is made between calcu- 
lations with 51 grid points and with 76 grid points for the initial velocity 

magnitude 0.175. The results show that the effect of the increment in number 
of grid points is almost negligible, and imply the convergence of the finite 
difference scheme. As for round-off error, calculations with double-precision 
accuracy were made at the wave magnitude 0.7 and the results are shown in the 
following table. After 243 calculation steps, the difference is still 0(10~5) 
and very small. Thus the finite difference scheme is very stable. 
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Table 3: Effect of round-off error; single precision (accuracy of 7 deci- 

mals) and double precision (accuracy of 16 decimals) . 


Time Steps 

Position 

(Double) 

Time 

Velocity 

Pressure 


0 



Initial Condition 




53 

0.488473 

(0.488490) 

1.05907 

(1.05909) 

0.544833 

(0.544851 

X 

X 

io -i 
10 ) 

-0.445474 

(-0.445476 

X 

X 

io -i 
10 ) 

103 

0.474543 

2.06016 

0.510547 

X 

10_ 2 

0.426178 

X 

10 ~! 


(0.474523) 

(2.06021) 

(0.511077 

X 

10- 2 ) 

(0.426278 

X 

10 ) 

143 

0.482243 

3.06948 

-0.643517 

X 

io ’i 

-0.397105 

X 

io l! 


(0.482259) 

(3.06952) 

(-0.643588 

X 

10 x ) 

(-0.397099 

X 

10 ) 

193 

0.474554 

4.03490 

0.997936 

X 

10 ~1 

0.213518 

X 

10-} 


(0.474553) 

(4.03493) 

(0.998079 

X 

10 ) 

(0.213533 

X 

10 x ) 

243 

0.473359 

5.03056 

-0.923763 

X 

10 ~! 

0.313444 

X 

1( T} 


(0.473345) 

(5.03058) 

(-0.923878 

X 

10 ) 

(0.313555 

X 

10 ) 


Case (b) . - For oscillations driven by wall amplification (eq. (27)), 
six of a values were chosen and the results are shown in Fig. 11-A and B. 

Table 4: Shock occurrences for various rates of wall amplification. 


Shock 


Time 

Position 

Magnitude 

0.20 

5.1 

0.14 

0.356 

0.15 

6.1 

0.92 

0.307 

0.10 

11.1 

0.74 

0.390 

0.05 

18.3 

0.003 

0.361 

0.03 

31.1 

0.07 

0.355 

0.01 

72.2 

0.08 

0.238 
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A shock is most likely to occur at the center, but now the boundary condi- 
tion v = - op is introduced at the outer surface. Thus, when the wave is mov- 
ing inward towards the center of the cylinder, it is at the same time being 
pulled back at the wall, decreasing the energy flow toward the center. Thus, 
until the strength of the wave gets big enough to overcome this pullback, a 
shock will not be generated. That is the reason why the shock is developed 
at such a high magnitude in this case in contrast to cases (a) and (c) below. 

When a is below 0.05, a weak shock is generated almost at the center of 
the cylinder and is different from the ones when a is above 0.1. Figs. 12 and 
13 show the wave patterns for a = 0.05 and a = 0.1 just before the shock is 
generated. Also Fig. 14 shows- the process of wave amplification for these two 
cases . 

In the case a = 0.01, the amplification of wave magnitude is very small, 
as can be seen in Figs. 11-A and B. But still a shock is generated. This 
implies that the power series converges up to a certain value of e, and this 
e is no bigger than 0.090 (magnitude of velocity = 0.2). 

Case (c) . - In this case, the value a = 0.05 is chosen and the amplifica- 
tion is stopped at various times but always when most of the kinetic energy 
is transferred into potential energy. Thus, all actions of stopping amplifi- 
cation occur at the same phase of the oscillation. 

Figs. 15 and 16 show the results. In Fig. 16 time is measured from the 
time when the amplification is stopped. Thus, t = 0 line in Fig. 16 is 
equivalent to the a = 0.05 line in Fig. 15. 

Compared to Case (a), this case has much less fluctuation of magnitude, 
as was anticipated. From Fig. 16, one can see that the largest magnitude at 
which a transverse wave canoscillate in the cylinder without growing into a 
shock is somewhwere around 0.13 (e = 0.058). This figure .is a little smaller 
than the result obtained for Case (a) (0.18). Thus the history of the oscill- 
ation must be an important factor in this problem. 

For all cases (a,b, and c) , it appears that there exists a permanent 
potential form for this resonant mode of transverse oscillation in the circu- 
lar cylinder for sufficiently small but finite values of e. 

For the purpose of checking whether the above results are sensitive to 
the value of a, one calculation was made stopping amplification around 
velocity magnitude 0.2 for the case of a = 0.01. The result is shown in 
Fig. 11-B. A shock is still generated in this case at the position R = 0.14. 

So we may conclude that the above results are quite independent of the value 
of a, as long as a is reasonably small. 
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CONCLUSION 


Numerical solutions have been obtained for the one-dimensional resonant 
mode (n = 0) of the nonlinear oscillation transverse to the axis of a circu- 
lar cylinder. 

Three types of calculations are performed, and each indicates the 
existence of a permanent potential form for this resonant mode (n = 0) . There 
is apparently a critical value of e for the power series to converge. The 
value is around 0.06 2 0.08. Above this value, the wave develops a shock 
pattern, but below it, the wave oscillates in the cylinder for an indefinite 
time, probably permanently, as suggested in Refs. 3 and 4. 
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APPENDIX A 


Calculation of Integral 


,D 


1 , , ,D 1 , 

— vdt and /_ — vdt 
r B r 



lines AD and BD. 


On the line AD, we have 


dr = t 1 + v] AD dt * 


Thus, the integration becomes, 


,D 1 , 1 

/. — vdt = 77 — -r 
A r [1+v] 


AD 


/"II r 
A r 


(A-l) 


Assuming linear velocity profile along the characteristic line AD, we can 
evaluate the integral. The result is, 


r D i , 1 

A r [T+v] 


AD 


r A ( W 

» V D - V + ' V A - 


In 0 } 
A 


* In the case when point A is at the center of the cylinder, i.e. r A = 0, 
we have 


v_ 

,D 1 D 

/. — vdt = r 1 . 

A r [1+v]^ 

In a similar way, we can integrate the equation on the line BD. 


* In the case when point D is at the center of the cylinder, i.e., r Q = 0, 
we have , 


V 

r D 1 , B 

/ — vdt = -jq r 

B r [1-vj 


BD 
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APPENDIX B 


The equation is 

<j) - (1 - <j> = 2d> 6 = .— <f) 

y tt v Y r y rr T r y rt r Y r 

subject to 4> (1 , t) = <j> r (0 , t ) = 0. 

We expand the velocity potential in a power series in amplitude 
parameter e : 

<j> = + e 2 4> ( ' 2 ' ) + (B-l) 


Then, the first order equation is 


(1) _ ,(D 


r tt 


rr 


i 

r 


= 0 


with the following boundary conditions: 


L 1} (l,t) = d)5 1) (0,t) = 


From the boundary conditions, 


= g(t) j Q (Br) 


where B is a constant such that (8) = 0, i.e., B = 3.83171. Also, g(t) 
will be cos (Bt) or sin (Bt) , so we can get a solution in the form, 


<1>^ = cos (Bt) J Q (3r) 


(B-2) 


The second order equation is 


(2) _ (2) _ 1 (2) = _ (1) (1) 

tt rr r T r Y r y rt 
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Introducing eq. (B-2) into the above equation yields 


(2) _ ,(2) _ 1 ,(2) „ o 3 T 2 


^tt' ' ^rr' " 7 * t ' = 3 J 1 (Br) Sin (2Bt) 


By the same technique of separation of variables, 


= 3f(r) sin (23t), 


(B-3) 


and, for the function f(r), 


f + - f + 4g 2 f = 
rr r r 


3 2 J 2 (Br) 


(B-4) 


with the boundary conditions: 


f r (0) = f r (l) = 0 


Using a Green' s-Funct ion method, the inhomogeous Bessel's equation (B-4) 
is solved: 


£ l - J 0 (2Sr) 

f R = j 0 ^ 2Br ^ ” Y (2 3) T 0 


The Wronskian is 


_2_ J(2g) 
■nr Y (2g) 


Then 

_ g Y (23) 

f = 2 {/ 0 x J q (2x) J 2 (x) dx y ~ ( 2 g ) ~ J Q ( 2 3r) 

- / Br x J q (2x) J 2 (x) dx Y Q (23r) 

- / B r x J 2 (x) Y q (2x) dx J Q (23r) } (B-5) 
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The integrations in eq. (B-5) are carried out by using the trapezoidal rule. 
Thus the solution of the original equation, to second order, is, 


<j> = e cos (gt) Jq (gr) + e 2 g sin (2gt) f(r) (B-6) 

where f(r) is given by Eq. (B-5). Then initial conditions (at t = 0) for u 
and v are 


u = 4> t (t = 0) = 2 (ge) 2 f(r) 

v = <j> (t = 0) = - (eg) J (gr) 
r 1 


(B-7) 


and the magnitude of the initial wave pattern is defined by 


Mag. = 0.581864 (eg). 


Also, the magnitude of wave over a certain cycle is defined by the largest 
amplitude of the wave during the period. 
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Figure 7. Shock occurrences for various initial wave magnitudes of case (a). Here, 
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Figure 8. Estimate of the initial wave magnitude at which oscillation 
is permanent, for case (a). 


24 


(°D/A) A 


0 

1 






27 


Figure 10. Convergence test for the wave of case (a) with initial velocity magnitude 0.175. 
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Figure 11-A. Wave amplifications and shock occurrences for case (b). Boundary condition 
at the wall is: v = -ap . 



TIME ( aot/Ro) 

Figure 11-B. Wave amplifications and shock occurrences for case (b). Boundary condition at the wall is: 
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Figure 12. Wave pattern just before'shock, for case (b) 
with amplification coefficient a = 0.1. 

Shock is generated at time =11.1 and 
position = 0.74. 
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Figure 13- Tlg^ aid position 
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Figure 14. Wave amplification process for case 
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various velocity magnitudes. R in bracket indicates non-dimensional 
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NASA-Langley, 1972 12 


F-7114 


Figure 16. Estimate of the wave magnitude at which oscillation is permanent, for case (c) 
Ordinate (Time = 0 line) corresponds to the wave growth line with the 
amplification coefficient o = 0.05 in Figure 15, i.e., time is measured from 
the time when the amplification is stopped. Circled points are transferred 
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